subroutine move_object(lid,sid,px,py,pz)
  use param
  use event_param
  use table_param,only:sim_obj
  use table_param,only:alloc_event_n_in_sim_obj
  use table_param,only:event_table
  use table_param,only:n_event_table
  use table_param,only:del_event_table
  use table_param,only:n_del_event_table
  use table_param,only:event_profile
  use table_param,only:n_event_profile
  use table_param,only:alloc_event_n_in_profile
  use table_param,only:spont_event_table
  use table_param,only:n_spont_event_table
  use box_param,only:vbox
  use box_param,only:vbox_map
  use box_param,only:cell_width
  use box_param,only:ch_width
  use label_data_param,only:obj_type
  use react_data_param,only:react_type
  use react_data_param,only:react_map
  use geom_data_param,only:geom_obj
  use geom_data_param,only:n_geom_obj
  use sink_data_param,only:sink_type
  use sink_data_param,only:sink_map
  use count_param,only:count_obj
  use profile_map_param,only:profile_map_rctn
  use profile_map_param,only:profile_map_sink
  use global_param,only:n_point_type
  implicit none

  interface
     function reallocInt(p,n)
       integer,pointer,dimension(:) :: reallocInt,p
       integer,intent(in) :: n
     end function reallocInt
  end interface
  
  integer,intent(in) :: lid,sid
  real(dbl),intent(in) :: px,py,pz

  integer :: i,j,k
  integer :: cur,bid,gid,nid,eid,pid,rid
  integer :: sid2,eid2,srank,ne,etype,nr,prank,event_id
  integer :: nx0,ny0,nz0
  integer :: n,nx,ny,nz,npos
  real(dbl) :: hx,hy,hz,dist2

  ! Assign the positions to sim_obj
  sim_obj(sid)%pos(1)=px
  sim_obj(sid)%pos(2)=py
  sim_obj(sid)%pos(3)=pz

  ! Remove it from the virtual box
  bid=sim_obj(sid)%vbox_index
  if (sim_obj(sid)%prev.eq.sid) then
     gid=sim_obj(sid)%next
     vbox(bid)%first=gid
     sim_obj(gid)%prev=gid
  else if (sim_obj(sid)%next.eq.sid) then
     gid=sim_obj(sid)%prev
     sim_obj(gid)%next=gid
  else
     sim_obj(sim_obj(sid)%prev)%next=sim_obj(sid)%next
     sim_obj(sim_obj(sid)%next)%prev=sim_obj(sid)%prev
  end if
  vbox(bid)%n_obj=vbox(bid)%n_obj-1

  ! Index for the new vbox
  nx=ceiling(px/cell_width)
  ny=ceiling(py/cell_width)
  nz=ceiling(pz/cell_width)

  bid=vbox_map(nx,ny,nz)
  sim_obj(sid)%vbox_index=bid
  sim_obj(sid)%next=0
  sim_obj(sid)%prev=0

  ! Count newly
  ! Old indices for the count box
  nx0=sim_obj(sid)%box(1)
  ny0=sim_obj(sid)%box(2)
  nz0=sim_obj(sid)%box(3)

  ! New indices for the count box
  nx=ceiling(px/ch_width(1))
  ny=ceiling(py/ch_width(2))
  nz=ceiling(pz/ch_width(3))

  do i=1,obj_type(lid)%n_group_belong
     gid=obj_type(lid)%group_list(i)
     nid=obj_type(lid)%config(count_obj(gid)%point_type)
     if (nx0.ne.nx) then
        if (count_obj(gid)%group_type.lt.3) then
           count_obj(gid)%nx(nx0)=count_obj(gid)%nx(nx0)-1
           count_obj(gid)%nx(nx)=count_obj(gid)%nx(nx)+1
        else
           count_obj(gid)%nx(nx0)=count_obj(gid)%nx(nx0)-nid
           count_obj(gid)%nx(nx)=count_obj(gid)%nx(nx)+nid
        end if
     end if
     if (ny0.ne.ny) then
        if (count_obj(gid)%group_type.lt.3) then
           count_obj(gid)%ny(ny0)=count_obj(gid)%ny(ny0)-1
           count_obj(gid)%ny(ny)=count_obj(gid)%ny(ny)+1           
        else
           count_obj(gid)%ny(ny0)=count_obj(gid)%ny(ny0)-nid
           count_obj(gid)%ny(ny)=count_obj(gid)%ny(ny)+nid
        end if
     end if
     if (nz0.ne.nz) then
        if (count_obj(gid)%group_type.lt.3) then
           count_obj(gid)%nz(nz0)=count_obj(gid)%nz(nz0)-1
           count_obj(gid)%nz(nz)=count_obj(gid)%nz(nz)+1           
        else
           count_obj(gid)%nz(nz0)=count_obj(gid)%nz(nz0)-nid
           count_obj(gid)%nz(nz)=count_obj(gid)%nz(nz)+nid
        end if
     end if
  end do

  sim_obj(sid)%box(1)=nx
  sim_obj(sid)%box(2)=ny
  sim_obj(sid)%box(3)=nz

  ! Remove only the associated reaction events from the event table
  nr=0
  do i=1,sim_obj(sid)%n_event
     eid=sim_obj(sid)%event_table_list(i)
     etype=event_table(eid)%event_type
     select case(etype)
     case(rctn_event)
        ! Adjust the event table for the other object
        if (sid.ne.event_table(eid)%obj1) then
           sid2=event_table(eid)%obj1
           srank=event_table(eid)%obj_seq1
           ne=sim_obj(sid2)%n_event
        else
           sid2=event_table(eid)%obj2
           srank=event_table(eid)%obj_seq2
           ne=sim_obj(sid2)%n_event
        end if
        eid2=sim_obj(sid2)%event_table_list(ne)
        sim_obj(sid2)%event_table_list(srank)=eid2
        ne=ne-1
        sim_obj(sid2)%n_event=ne
        if (sid2.eq.event_table(eid2)%obj1) then
           event_table(eid2)%obj_seq1=srank
        else
           event_table(eid2)%obj_seq2=srank
        end if

        ! Remove the event table
        n_del_event_table=n_del_event_table+1
        del_event_table(n_del_event_table)=eid
        event_table(eid)%removed=.true.
        n_event_table=n_event_table-1

        ! To update the n_event in sim_obj
        nr=nr+1
        
        ! Update the event profile
        pid=event_table(eid)%profile_table_index
        prank=event_table(eid)%profile_table_seq
        ne=event_profile(pid)%n_event
        eid2=event_profile(pid)%event_list(ne)
        event_profile(pid)%event_list(prank)=eid2
        event_table(eid2)%profile_table_seq=prank
        ne=ne-1
        event_profile(pid)%n_event=ne
     case(sink_event)
        ! Remove the event from the table
        n_del_event_table=n_del_event_table+1
        del_event_table(n_del_event_table)=eid
        event_table(eid)%removed=.true.
        n_event_table=n_event_table-1

        ! To update the n_event in sim_obj
        nr=nr+1

        ! Update the event profile
        pid=event_table(eid)%profile_table_index
        prank=event_table(eid)%profile_table_seq
        ne=event_profile(pid)%n_event
        eid2=event_profile(pid)%event_list(ne)
        event_profile(pid)%event_list(prank)=eid2
        event_table(eid2)%profile_table_seq=prank
        ne=ne-1
        event_profile(pid)%n_event=ne
     end select
  end do

  ! Adjust the number of events belonging to the current sim_obj
  sim_obj(sid)%n_event=sim_obj(sid)%n_event-nr

  ! Save the index of event_table_list in sim_obj
  event_id=sim_obj(sid)%n_event

  ! Determine the reaction events newly
  do i=1,vbox(bid)%n_neighbor

     ! Sweeping over the neighbor cells
     nid=vbox(bid)%neighbor_list(i)
     cur=vbox(nid)%first

     ! Loop over the objects in the cell
     do j=1,vbox(nid)%n_obj
        
        ! Check if the reaction is defined in the data sheet.
        rid=react_map(lid,sim_obj(cur)%object_type)
        if (rid.eq.0) go to 100

        ! Calculate the distance between objects
        hx=sim_obj(cur)%pos(1)+vbox(bid)%dx(i)-px
        hy=sim_obj(cur)%pos(2)+vbox(bid)%dy(i)-py
        hz=sim_obj(cur)%pos(3)+vbox(bid)%dz(i)-pz
        dist2=hx*hx+hy*hy+hz*hz

        ! Loop over the possible reaction types
        do k=1,react_type(rid)%n_data
           
           ! Compare the distance and capture radius
           if (dist2.gt.react_type(rid)%react_rd2(k)) cycle

           ! Determine the index of the event table
           n_event_table=n_event_table+1
           call alloc_event_table(n_event_table,eid)

           ! Set the object 1
           event_id=event_id+1
           sim_obj(sid)%n_event=event_id
           if (event_id.gt.sim_obj(sid)%n_alloc_event) then
              sim_obj(sid)%n_alloc_event=&
                   sim_obj(sid)%n_alloc_event+&
                   alloc_event_n_in_sim_obj
              sim_obj(sid)%event_table_list=>&
                   reallocInt(sim_obj(sid)%event_table_list,&
                   sim_obj(sid)%n_alloc_event)
           end if
           sim_obj(sid)%event_table_list(event_id)=eid
           
           ! Set the object 2
           sim_obj(cur)%n_event=sim_obj(cur)%n_event+1
           if (sim_obj(cur)%n_event.gt.&
                sim_obj(cur)%n_alloc_event) then
              sim_obj(cur)%n_alloc_event=&
                   sim_obj(cur)%n_alloc_event+&
                   alloc_event_n_in_sim_obj
              sim_obj(cur)%event_table_list=>&
                   reallocInt(sim_obj(cur)%event_table_list,&
                   sim_obj(cur)%n_alloc_event)
           end if
           sim_obj(cur)%event_table_list(sim_obj(cur)%n_event)=eid

           ! Set the event table
           event_table(eid)%obj1=sid
           event_table(eid)%obj_seq1=event_id
           event_table(eid)%obj2=cur
           event_table(eid)%obj_seq2=sim_obj(cur)%n_event
           event_table(eid)%event_type=rctn_event
           event_table(eid)%sqnc=k
           event_table(eid)%react_index=rid
           event_table(eid)%removed=.false.
           event_table(eid)%is_spontaneous=.false.

           ! Set the event profile
           if (react_type(rid)%react_const_rate(k).eq.0) then
              ! if the event is spontaneous,
              ! add to the spont_event_table.
              event_table(eid)%is_spontaneous=.true.
              n_spont_event_table=n_spont_event_table+1
              spont_event_table(n_spont_event_table)=eid
              event_table(eid)%spont_table_index=&
                   n_spont_event_table
           else
              ! if the event is probabilistic,
              ! add to the event profile
              if (profile_map_rctn(rid,k).eq.0) then
                 ! Create a new profile
                 n_event_profile=n_event_profile+1
                 call alloc_profile_table(n_event_profile,pid)
                 profile_map_rctn(rid,k)=pid
                 event_profile(pid)%event_type=rctn_event
                 event_profile(pid)%rate_const=&
                      react_type(rid)%react_const_rate(k)
                 event_profile(pid)%index=rid
                 event_profile(pid)%sqnc=k
                 event_profile(pid)%n_event=1
                 if (event_profile(pid)%n_alloc_event.eq.0) then
                    event_profile(pid)%n_alloc_event=&
                         alloc_event_n_in_profile
                    event_profile(pid)%event_list=>&
                         reallocInt(event_profile(pid)%event_list,&
                         event_profile(pid)%n_alloc_event)
                 end if
                 event_profile(pid)%event_list(1)=eid
                 event_table(eid)%profile_table_index=pid
                 event_table(eid)%profile_table_seq=1
              else
                 ! Add to the existing profile
                 pid=profile_map_rctn(rid,k)
                 n=event_profile(pid)%n_event
                 n=n+1
                 event_profile(pid)%n_event=n
                 if (n.gt.event_profile(pid)%n_alloc_event) then
                    event_profile(pid)%n_alloc_event=&
                         event_profile(pid)%n_alloc_event+&
                         alloc_event_n_in_profile
                    event_profile(pid)%event_list=>&
                         reallocInt(event_profile(pid)%event_list,&
                         event_profile(pid)%n_alloc_event)
                 end if
                 event_profile(pid)%event_list(n)=eid
                 event_table(eid)%profile_table_index=pid
                 event_table(eid)%profile_table_seq=n
              end if
           end if
        end do
100     cur=sim_obj(cur)%next
     end do
  end do

  ! Sink with the geometric object
  do i=1,n_geom_obj
     
     ! Check if the reaction is defined in the data sheet
     gid=sink_map(i,lid)
     if (gid.eq.0) cycle

     ! Calculate the distance between the geom_obj and the object
     select case(geom_obj(i)%geom_type)
     case(1)
        ! point
        hx=px-geom_obj(i)%x1
        hy=py-geom_obj(i)%y1
        hz=pz-geom_obj(i)%z1
        dist2=hx*hx+hy*hy+hz*hz
     case(2)
        ! straight line
        npos=((geom_obj(i)%x1-px)*(geom_obj(i)%x1-geom_obj(i)%x2)&
             +(geom_obj(i)%y1-py)*(geom_obj(i)%y1-geom_obj(i)%y2)&
             +(geom_obj(i)%z1-pz)*(geom_obj(i)%z1-geom_obj(i)%z2))/&
             ((geom_obj(i)%x1-geom_obj(i)%x2)**2&
             +(geom_obj(i)%y1-geom_obj(i)%y2)**2&
             +(geom_obj(i)%z1-geom_obj(i)%z2)**2)
        if (npos.lt.0) then
           hx=geom_obj(i)%x1
           hy=geom_obj(i)%y1
           hz=geom_obj(i)%z1
        else if (npos.gt.1) then
           hx=geom_obj(i)%x2
           hy=geom_obj(i)%y2
           hz=geom_obj(i)%z2
        else
           hx=npos*(geom_obj(i)%x2-geom_obj(i)%x1)+geom_obj(i)%x1
           hy=npos*(geom_obj(i)%y2-geom_obj(i)%y1)+geom_obj(i)%y1
           hz=npos*(geom_obj(i)%z2-geom_obj(i)%z1)+geom_obj(i)%z1
        end if
        dist2=(px-hx)**2+(py-hy)**2+(pz-hz)**2
     case(3)
        dist2=(px-geom_obj(i)%x1)**2
     case(4)
        dist2=(px-geom_obj(i)%x1)**2
     case(5)
        dist2=(py-geom_obj(i)%y1)**2
     case(6)
        dist2=(py-geom_obj(i)%y1)**2
     case(7)
        dist2=(pz-geom_obj(i)%z1)**2
     case(8)
        dist2=(pz-geom_obj(i)%z1)**2
     end select

     ! loop over the each type of sink reaction
     do j=1,sink_type(gid)%n_data

        ! Compare the distance and capture radius
        if (dist2.gt.sink_type(gid)%react_rd2(j)) cycle

        ! Determine the index of the event table
        n_event_table=n_event_table+1
        call alloc_event_table(n_event_table,eid)

        ! Set the object
        event_id=event_id+1
        sim_obj(sid)%n_event=event_id
        if (event_id.gt.sim_obj(sid)%n_alloc_event) then
           sim_obj(sid)%n_alloc_event=&
                sim_obj(sid)%n_alloc_event+&
                alloc_event_n_in_sim_obj
           sim_obj(sid)%event_table_list=>&
                reallocInt(sim_obj(sid)%event_table_list,&
                sim_obj(sid)%n_alloc_event)
        end if
        sim_obj(sid)%event_table_list(event_id)=eid

        ! Set the event table
        event_table(eid)%obj1=sid
        event_table(eid)%obj_seq1=event_id
        event_table(eid)%obj2=i
        event_table(eid)%obj_seq2=j
        event_table(eid)%event_type=sink_event
        event_table(eid)%react_index=gid
        event_table(eid)%removed=.false.
        event_table(eid)%is_spontaneous=.false.

        ! Set the event profile
        if (sink_type(gid)%react_const_rate(j).eq.0) then
           ! if the event is spontaneous,
           ! add to the spont_event_table
           event_table(eid)%is_spontaneous=.true.
           n_spont_event_table=n_spont_event_table+1
           spont_event_table(n_spont_event_table)=eid
           event_table(eid)%spont_table_index=&
                n_spont_event_table
        else
           ! if the event is probabilistic,
           ! add to the event_profile
           if (profile_map_sink(gid,j).eq.0) then
              ! Create a new profile
              n_event_profile=n_event_profile+1
              call alloc_profile_table(n_event_profile,pid)
              profile_map_sink(gid,j)=pid
              event_profile(pid)%event_type=sink_event
              event_profile(pid)%rate_const=&
                   sink_type(gid)%react_const_rate(j)
              event_profile(pid)%sqnc=j
              event_profile(pid)%n_event=1
              if (event_profile(pid)%n_alloc_event.eq.0) then
                 event_profile(pid)%n_alloc_event=&
                      alloc_event_n_in_profile
                 event_profile(pid)%event_list=>&
                      reallocInt(event_profile(pid)%event_list,&
                      event_profile(pid)%n_alloc_event)
              end if
              event_profile(pid)%event_list(1)=eid
              event_table(eid)%profile_table_index=pid
              event_table(eid)%profile_table_seq=1
           else
              ! Add to the existing profile
              pid=profile_map_sink(gid,j)
              n=event_profile(pid)%n_event
              n=n+1
              event_profile(pid)%n_event=n
              if (n.gt.event_profile(pid)%n_alloc_event) then
                 event_profile(pid)%n_alloc_event=&
                      event_profile(pid)%n_alloc_event+&
                      alloc_event_n_in_profile
                 event_profile(pid)%event_list=>&
                      reallocInt(event_profile(pid)%event_list,&
                      event_profile(pid)%n_alloc_event)
              end if
              event_profile(pid)%event_list(n)=eid
              event_table(eid)%profile_table_index=pid
              event_table(eid)%profile_table_seq=n
           end if
        end if
     end do
  end do

  ! Add to the corresponding virtual box
  if (vbox(bid)%n_obj.eq.0) then
     vbox(bid)%first=sid
     sim_obj(sid)%next=sid
     sim_obj(sid)%prev=sid
  else
     sim_obj(sid)%next=vbox(bid)%first
     sim_obj(vbox(bid)%first)%prev=sid
     vbox(bid)%first=sid
     sim_obj(sid)%prev=sid
  end if
  vbox(bid)%n_obj=vbox(bid)%n_obj+1

  return
  
end subroutine move_object
